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For some years Lanczos moments methods have been combined with large-scale shell-model cal- 
culations in evaluations of the spectral distributions of certain operators. This technique is of great 
value because the alternative, a state-by-state summation over final states, is generally not feasible. 
The most celebrated application is to the Gamow- Teller operator, which governs j3 decay and neu- 
trino reactions in the allowed limit. The Lanczos procedure determines the nuclear response along 
a line q=0 in the {to, q) plane, where u> and q are the energy and three-momentum transfered to 
the nucleus. However, generalizing such treatments from the allowed limit to general electroweak 
response functions at arbitrary momentum transfers seems considerably more difficult: the response 
function must be determined over the entire (u),q) plane for an operator 0(q) that is not fixed, but 
depends explicitly on q. Such operators arise in any semileptonic process where the momentum 
transfer is comparable to (or larger than) the inverse nuclear size. Here we show, for Slater deter- 
minants built on harmonic oscillator basis functions, that the nuclear response for any multipole 
operator O(q) can be determined efficiently over the full response plane by a generalization of the 
standard Lanczos moments method. We describe the Piecewise Moments Method and thoroughly 
explore its convergence properties for the test case of electromagnetic responses in a full sd-shell 
calculation of 28 Si. We discuss possible extensions to a variety of electroweak processes, including 
charged- and neutral-current neutrino scattering. 



I. INTRODUCTION AND MOTIVATION 

The shell model [| is an important tool for modeling 
the ground states and low-lying excitations of nuclei, as 
well as their electromagnetic and weak-interaction prop- 
erties. Calculations are performed by diagonalizing an 
effective interaction H within a low-momentum model 
space typically consisting of all Slater determinants that 
can be constructed within one or perhaps a few princi- 
pal shells. H corrects for the absence of the remaining 
shells and of the high-momentum correlations that reside 
primarily in those excluded shells. Typically H is deter- 
mined empirically 0] , though recently considerable work 
has been invested in various effective-theory approaches 
to derive H directly from the underlying bare interaction 
0. 

With increasing CPU speeds and the advent of paral- 
lel computing, shell model practitioners have been able 
to tackle problems in very large bases - ranging in some 
cases to Hamiltonians of dimension ~ 10 s — 10 9 pj. Be- 
cause direct diagonalizations in such spaces are impossi- 
ble, iterative methods are used to determine extremum 
eigenvalues and their eigenfunctions. The Lanczos al- 
gorithm is the most commonly used method. It is 
based on a recursive mapping of the full Hamiltonian 
into a much smaller, tridiagonal matrix that preserves 
certain exact information on the lowest moments of the 
full Hamiltonian. As described in the next section, ex- 
tremum eigenvalues of the Lanczos matrix converge to 
those of the full Hamiltonian because of the algorithm's 
moments properties, thus making the algorithm useful 



as a diagonalization tool. However, its real power is con- 
nected with another property, a stable method for solving 
the classical moments problem - determining the sim- 
plest discrete distribution characterized by those same 
lowest moments. 

One can argue that the Lanczos algorithm adresses 
both of the most common challenges encountered in nu- 
clear structure physics - detailed information about spe- 
cific low-lying states and global information on spectral 
moments important to inclusive properties, such as re- 
sponse functions, polarizabilities, and Green's functions. 

The moments method, described in more detail in the 
next section, can be applied in a straight-forward way to 
operators like the Gamow- Teller (GT) operator, 
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This operator is independent of the three-momentum 
transfer q. GT strength distributions (and their neutral- 
current analogs) govern low-energy neutrino reactions 
important to solar and supernova physics. The combi- 
nation of increasingly sophisticated shell-model ground- 
state wave functions and Lanczos methods for strength 
distributions has had great impact. Results include much 
better agreement between theory and experiment (e.g., 
(p, n) mappings of GT distributions) and new iron-group 
weak interaction rates [j| that have influenced the col- 
lapse and explosion physics of Type II supernovae 0, Q • 
However, the restriction to allowed operators is very 
limiting. For example, in supernova physics, approxi- 
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mately half of the strength for heavy-flavor neutral cur- 
rent scattering is carried by momentum-dependent oper- 
ators JJ. Because no efficient moments method is avail- 
able for these more complex operators, and because state- 
by-state summations would be prohibitively difficult, the- 
orists have not been able to make use of state-of-the-art 
shell model wave functions in evaluating inclusive forbid- 
den response functions. Instead less sophisticated meth- 
ods - simplified shell model treatments, QRPA, or even 
schematic approaches like the Goldhaber- Teller model - 
have been substituted, as the resulting smaller Hilbcrt 
spaces do allow state-by-state summations. The difficulty 
in designing a moments method for a multipole operator 
0(q) that varies with q is clear: the standard moments 
method would require one to repeat the calculation many 
times over a grid in q, as 0(q) is a fixed operator only 
along g=constant trajectories. Furthermore, as will be- 
come clearer below, a rather dense grid in q would be 
needed to interpret an experiment that maps out some 
nontrivial trajectory in the (w, q) response plane, as the 
discreteness of the moments distributions at each q makes 
it complicated to extrapolate in q. 

Interest in general electroweak response functions is 
not limited to supernova physics, clearly. An electron 
scattering experiment, for example, generally maps out 
some area within the space-like half of the response plane, 
depending on the range of spectrometer angles and elec- 
tron energy loss explored. 

In this paper we show that there is an efficient mo- 
ments technique for constructing the entire shell-model 
inclusive response surface - the response S(cu, q) as a con- 
tinuous function of both q and u, as well as of the oscilla- 
tor parameter b - provided the underlying single-particle 
basis wave functions are taken to be harmonic oscillators. 
The method involves a small number of Lanczos calcu- 
lations (at most three in the 28 Si test case we use here). 
That is, the full response surface can be reconstructed 
with little more effort than required in the familiar GT 
case, where results are obtained only for a line q=0 in the 
response plane. The information that must be extracted 
from the many-body calculations to perform this recon- 
struction are the elements of the tridiagonal Lanczos ma- 
trices and certain dot products between Lanczos vectors. 
The method can be viewed as a nearly perfect "numeri- 
cal effective theory," as it extracts from potentially quite 
complicated microscopic calculations precisely that infor- 
mation necessary to reconstruct the response surface up 
to a specified resolution in energy This has implications 
for problems like modeling supernova explosions, where 
one issue in the use of sophisticated nuclear physics re- 
sults is defining a practical scheme for making the infor- 
mation available "on line" within the explosion codes. 

The outline of this paper is as follows: In the next 
section we review in more detail the Lanczos method, 
including its conventional uses in constructing response 
functions and Green's functions. In the third section we 
discuss the problem at hand, the construction of response 
functions for electroweak nuclear operators at arbitrary 



(w,q). We discuss properties of harmonic oscillator ma- 
trix elements of these operators which could lead to ef- 
ficient algorithms for constructing the full response sur- 
face. We describe numerical limitations to some possible 
approaches. In the fourth section we describe the Piece- 
wise Moments Method (PMM), which we formulate first 
for a Lorentzian resolution function in energy, but then 
generalize for any choice. In the fifth section we explore 
convergence properties of the PMM, using as a test case 
the electromagnetic response functions for a full sd-shell 
calculation of 28 Si. Convergence properties of the method 
are explored and various numerical comparisons are made 
with "exact" results. A series of results for response sur- 
faces are presented. In the concluding section we discuss 
opportunities for applying the method, including super- 
nova physics and the neutrino reactions at energies below 
1 GeV (e.g., LSND or atmospheric neutrinos). We sug- 
gest extensions of the work - to check unitarity and to 
project spurious states - that might be important in fu- 
ture applications of the PMM. 

II. LANCZOS ALGORITHM PRELIMINARIES 

The Lanczos method is based on a mapping of a Hamil- 
tonian H of dimension N to tridiagonal form by a recur- 
sive definition of a new orthonormal basis \vi), called the 
Lanczos vectors. One begins with an arbitrary unit vec- 
tor \vi) and performs the successive operations to define 
the \vt), 

H\v 1 )= «l|«l)+AI«2> 

H\v 2 ) = Pi\vi) + a 2 \v 2 ) + /3 2 \v 3 ) 

H\v 3 ) = (3 2 \v 2 ) + a 3 \v 3 ) + /3 3 \v 4 ) 

H\vi)= p i - 1 \v i - 1 )+a i \v i )+l3 i \v i+1 ). (1) 

where in the first step ai|i>i) is the projection of H\vi) 
onto \vi), while the remaining orthogonal portion defines 
a new unit vector and amplitude Pi\v 2 }- In the third 
step we see the tridiagonal form emerge, as the construc- 
tion demands (v 3 \H\vi) = 0. The hermiticity of H has 
been used above. In practice this algorithm, when used 
to determine extremum eigenvalues and eigenvectors, is 
usually executed with a reorthogonalization step in each 
iteration, to guarantee that the new vector \vi) is or- 
thogonal to all previous vectors. While this step is not 
required mathematically, it is nevertheless important nu- 
merically, as roundoff error can grow with successive it- 
erations [3 , eventually leading to a Lanczos matrix that 
contains multiple copies of subspaces of H. 

Clearly, if the procedure were continued N steps, one 
would find (3n =0, as exhaustion of the Hubert space 
terminates the construction. The effect of the Lanc- 
zos procedure would be a unitary transformation to a 
new basis in which the Hamiltonian is tridiagonal, but 
still of dimension N . However, in useful applications 
N is very large, and the Lanczos construction is trun- 
cated by choice after n iterations by setting /3 n = 0, 
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n <C N. This yields the truncated tridiagonal Lanczos 
matrix L(n, 



L(n, \vx}) 
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The notation emphasizes that L is entirely determined 
by the number of iterations n and by the choice of the 
starting vector \v\). 

The powerful property of this mapping of H onto a 
much smaller subspace is that it extracts specific exact 
information from the full Hamiltonian 



H(N) 



{(vi\H x \ Vl ) 
{an, ■ ■ 



L(n, \vi)) 
X = l,...,2n- 

i Pi j ••*) ftn — 



1} 
} 



(3) 



That is, the first 2n— 1 moments of H for the starting vec- 
tor \vi) determine L(n, \vi}), and conversely. In a very 
real sense, the Lanczos algorithm can be considered a nu- 
merical effective theory. It systematically extracts from 
the large matrix the long-wavelength information de- 
scribing the distribution of \vx) over the eigenspectrum of 
H, while leaving behind the high frequency information 
important to the detailed structure of this distribution, 
but not to any of its broad features. More precisely, if we 
denote by {(tp Ei , Ei),i = 1, . . . , N} the exact eigenener- 
gies/functions of H and by {(tp E ,,Ei) 7 i = 1, . . . ,n} the 
eigenenergies/functions of L(n, |i>i)), then 

N 

i=l 

n 

i=X 
n 

= ElV i s i (l)| 2 ^A=l,... ) 2n-l (4) 

i=l 

That is, the eigenvalues and eigenvectors of the Lanc- 
zos matrix L(n, \v±}) determines a set of n points Ei and 
weights that solves the classical moments problem - a 
discrete distribution whose moments reproduce those of 
the full matrix H . The weights are simply the squares 
of the first components of the respective Lanczos eigen- 
vectors. As Whitehead has emphased | lCX ll l| . the speed 
and numerical stability of this classical moments solution 
is a very special property of the Lanczos algorithm. 

The usual application of the Lanczos algorithm is in de- 
termining cxtrcmum eigenvalues, especially the ground 
state and first few excited states. After each iteration 
L(n, \vi)) can be diagonalized by standard techniques, 
determining eigenvalues and eigenfunctions. Because ex- 
tremum eigenvalues are heavily weighted in H n when n 



is large and often separated by gaps from the bulk of 
spectrum, these eigenvalues of L must quickly converge 
to the true eigenvalues of H. This convergence can be 
monitored numerically as the algorithm is executed. In 
typical shell model applications the ground state often 
converges in about 50 iterations, while the lowest 10 or 
so eigenvalues and eigenfunctions may be resolved in 200 
iterations 

But perhaps the most elegant applications of the Lanc- 
zos algorithm are in distribution functions or Green's 
functions, inclusive quantities that depend most directly 
on spectral moments, the long- wavelength information 
extracted from H. Consider the response function S(u>) 



N 



(5) 



where the sum extends over a complete set of states i of 
the full Hamiltonian H . Assume for the moment that 
O is a fixed operator, like the Gamow- Teller operator, 
independent of the three-momentum transfer q to the 
nucleus. We define a unit vector \vi) by 



0\g.s. 



CVi 



(0) 



where c is the overall strength (norm). Taking \vi) as the 
starting vector and competing n Lanczos iterations, one 
can form the distribution 



S n (uj) = \c\ 2 yj(o;-E i )\^ i (l)\ 



(7) 



i=i 



If one weights this distribution with oj x and integrates 
over u, one sees from Eq. |ljthat S n {ui) reproduces the 
lowest In — 1 integrated moments of the exact spectral 
distribution S(lu). 

Now, as experiments are done with finite resolution, 
spectral variations occurring at energy scales below that 
resolution are irrelevant. Furthermore, the discrete spec- 
trum of the shell model is itself an artifact of the use 
of a finite Hilbert space. It represents resonances in 
the continuum by discrete doorway states, with the den- 
sity of such states increasing as the shell-model space 
is expanded, to better represent the continuum. Thus 
one quickly recognizes that the high-frequency informa- 
tion missing from Eq. \7\ may be irrelevant when com- 
paring with experiment. It is customary to replace the 
^-function in the Lanczos strength distribution with a 
resolution function, choosing a width parameter a ap- 
propriate to the experiment in question 



S n (uj, v) = c 2 J2 ~ A, ^)|^(i)| 5 



(8) 



where R is normalized to 1. Common choices are 
Lorentzians and Gaussians, e.g., 

a 1 



Rl(u - Ei,a) 
R g (uj - Ei,a) 



7T ( w - E,y + a 2 

' :exp(-(o;-^) 2 /2a) (9) 
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Because such smoothing makes high-frequency informa- 
tion irrelevant, the Lanczos approximation S n (uj,a) con- 
verges to the exact smoothed distribution S(u>,cr) for n 
sufficiently large. Qualitatively, convergence is achieved 
when the typical separation in energy of neighboring 
Lanczos eigenvectors becomes substantially smaller than 
the chosen resolution a. For a ~ 0.25 MeV, a value 
typical of (p,n) mappings of GT strength, for instance, 
this may occur at n ~ 200. Thus the Lanczos algorithm 
provides an exact method for determining appropriately 
smoothed strength functions in very complex shell- model 
spaces. 

A second spectral application is to Green's functions 
(l2^ . which arise in nuclear calculations of polarizabili- 
ties, in virtual processes like double [3 decay, and in a 
variety of many-body applications, such as effective in- 
teractions and operators. The Green's function acting on 
a normalized vector \v\) 

G(,)h)~K) (10) 
oj — H 

can be approximated after n Lanczos iterations as 

G n (w)\vi) = Pi(w)|ui) + <72(w)|u 2 )H \-g„(u>)\v n ) (11) 

where the coefficients gi(u>) are finite continued fractions 
formed from the entries in the tridiagonal matrix. For 
example 

91 M = " J2 

u-oti 

'■■ (12) 

With each additional iteration, one additional Lanczos 
vector is added to the expansion, and each continued 
fraction increases in rank by one through the addition of 
a new a n +i and (3 n . As most Green's function applica- 
tions involve convolutions with relatively smooth opera- 
tors, often G n (w)\v\) becomes numerically equivalent to 
G{ui)\v\) after a few Lanczos iterations (~ 20) |l3[ . 

An important consequence of Eq. ^2 is that, once the 
Lanczos calculation is completed, the Green's function is 
known as a function of w. This will be important in the 
applications we discuss later. 



III. ELECTROWEAK RESPONSE FUNCTIONS 
AT ARBITRARY q 2 

The discussion of the previous section addressed the 
special case of a fixed operator, like the GT operator, 
which governs the weak nuclear response along the <7=0 
line in the (u>,q) response plane. However many elec- 
troweak processes of interest - intermediate-energy elec- 
tron or neutrino scattering, muon capture, etc. - involve 



appreciable three-momentum transfers (and the associ- 
ated excitation of radial modes in the nucleus). That is, 
the relevant response function is 

N 

S(u,q) = ^|<Vx|0(<?)|<?. S .)| 2 <K W - Ei) (13) 
t=i 

where 0(q) is (an assumed one-body) electroweak oper- 
ator that depends explicitly on q. If one naively applies 
the formalism of the preceding section, a new calculation 
would be needed for each desired q, because the opera- 
tor evolves with q. This would require tediously stepping 
over a grid of fixed qs, computing a Lanczos calculation 
for each value, to map the full surface above the response 
plane. 

Here we discuss procedures for evaluating S(u,q) very 
efficiently as a function of q (and u>) over the entire re- 
sponse plane, at the cost of only a few Lanczos calcula- 
tions. The approach depends on the assumption that the 
shell-model basis of Slater determinants has been formed 
from harmonic oscillator single-particle wave functions. 
This choice allows one to exploit attractive properties of 
the matrix elements of 0{q) between such wave functions. 

While we will delay details of the test application 
(electromagnetic response functions for 28 Si) to the next 
section, here we sketch the basic idea. One can write 
0(q)\g.s.) in second quantization 

^2(a\0(q)\l3)aia^g.s.). (14) 

a, 13 

where a and (3 represent a complete set of single-particle 
quantum numbers. For the choice of harmonic oscilla- 
tors, matrix elements of the standard charge, longitudi- 
nal, transverse electric, and transverse magnetic multi- 
poles can be evaluated in closed form, leading to 0, 

(a\Oj(q)\(3) = y V- K V*e-Vp a < 3 (y). (15) 

Here we have denoted the multipolarity of the operator 
by J, K=2(l) for normal(abnormal) parity operators, and 
y = (qb/2) 2 with b the oscillator parameter. The crucial 
point is that p(y) is a finite polynomial in y or q 2 . In the 
28 Si test case, the most complicated operator that arises 
has only three nonzero terms in p(y). 

We first go through a schematic argument to show how 
this y dependence might be exploited. Denoting the order 
of the polynomial p by m, it follows that 

0(q)\g.s.) = 

y {J-K)/2 e -y [( , o ^ + Ciy |^ + . . . + Cmy ™\ vT )} 

^y(J-K)/2 e -y ciy)lMy)) (16) 

using a notation analogous to Eq. |SJ with the strength Cj 
chosen to make a unit vector. For parity-conserving 
interactions and standard phase conventions, all quanti- 
ties can be taken as real, with the cs nonnegative. The 
\v{), of course, are not orthonormal. Similarly c(y) and 
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1 171 (y)) can be viewed as a y-dependent strength and unit 
vector. It follows that 

N 

S(u>,q) = y J - K e- 2 y\c{y)\ 2 Y^\^EM{y))\ 2 5{uj - E t ) 

i=i 

(17) 

where 

m 

\c(y)\ 2 \(M\Mv))\ 2 = E ^^y'+'Ki^x^Ki), 

i,fc=o 

so that the response function has a similar polynomial 
form. It also follows that moments of S(u>) have the 
form 




S^duj = y J - K e- 2 y\c(y)\ 2 Y,\^E,\vi(y)}\ 2 E*. 

i=l 



(19) 

These last two results simply state that if one had a com- 
plete set of N eigenvalues and eigenfunctions, each con- 
tributing transition probability would have a simple, an- 
alytical behavior in y. 

Of course, these results are only of academic interest: 
as we are assuming N is prohibitively large, a complete 
diagonalization is impossible. This leaves a much more 
interesting question: can we find an analog of Eq.0or|Sl 
an efficient Lanczos representation of S(u>,q), that also 
exploits the polynomial behavior of the response in yl If 
so, it would appear to be a practical way to construct the 
response over the entire (w, q) plane. 

We have explored several of the possibilities, uncover- 
ing some of the numerical pitfalls. Even the less success- 
ful methods are interesting conceptually, so we describe 
the approaches qualitatively below, reserving details for 
the Appendix. Finally, we describe and test the Piecewise 
Moments Method, demonstrating that it is effectively ex- 
act and stable for an arbitrary number of iterations. 

The Naive Moments Method: The most straight- 
forward approach to the above problem exploits the 
equivalence between the tridiagonal matrix L(n, \v\(y))) 
and the moments (vi(y)\H x \vi(y)) , A = 1, . . . , 2n — 1. If 
we had a method to determine L(n, \vi(y))) as a contin- 
uous function of y, Eq.[7|could be then be used at any y. 
The result would be a discrete distribution along any line 
j/=constant, that would evolve smoothly as y is changed. 

The moments equivalence shows this is possible as 

1 m 

{vi(v)\h x Mv)) = J2 n ij< c iy i+j > (2°) 

where f2 A - = (v\\H x \v-[) are the "mixed moments" of H. 
Clearly, by operating successively with H on each \v\) 
n times, the ft}, can be evaluated. This then defines 
the moments at any y, and thus in principle the exact 
L(n,\vi(y))). This is then an exact moments Lanczos 
description of the response surface: the distribution that 



results from diagonalizing L(n, \vi(y))) will correctly de- 
scribe the moments (vi(y)\H x \vi(y)} , A = 1, . . . , 2n — 1. 
We call this the Naive Moments Method (NMM). 

The catch is the "in principle" part: the inversion from 
moments to L(n, \vi(y)}) is equivalent to the classical mo- 
ments problem: finding a discrete distribution from its 
moments. While specific formulae for this inversion are 
given in the Appendix, the inversion becomes increas- 
ingly unstable with increasing n. Even with calculations 
done in 64-bit precision, the NMM can fail in fewer than 
10 iterations. 

The Legendre Polynomial Moments Method: 

As discussed in the Appendix, the rapid loss of precision 
with increasing n in the NMM - more precisely, in the 
inversion to determine the Lanczos matrix and thus the 
distribution - can be traced to the dominance of the ex- 
tremum eigenvalues in high-order moments like (iJ 2 " _1 ). 
This suggests reformulating the NMM in such a way that 
the basic physics is preserved - the simple polynomial de- 
pendence of moments on y and the use of this dependence 
to define moments for all y - while using combinations of 
moments that are better behaved near the extremums. 

A possible choice to replace {1, if, H 2 , ....} are the Leg- 
endre polynomials in H, {P (H), Pi(H), P 2 (H), ...}, with 
the energy range between the lowest and highest eigen- 
values mapped on to [-1,1], the usual range. These poly- 
nomials have the attractive feature that they achieve a 
fixed magnitude of 1 at the boundaries of the range. The 
specific algorithm we constructed is described in the Ap- 
pendix. While equivalent to the NMM mathematically, 
the recurrence relations for determining the a, , from 
the Legendre polynomial moments indeed proved to be 
more stable. In some applications 80 iterations could be 
performed with little loss of precision. However, as dis- 
cussed in the Appendix, it sometimes fails more quickly. 

Other orthogonal polynomials in H could be used. An 
interesting question is the possibility that some choice 
might further improve the stability. 

Legendre Coefficients Method: Rather than car- 
rying through the procedure of constructing the Lanczos 
matrix, diagonalizing and then constructing the strength 
distribution, it is possible to find the a truncated expan- 
sion of the strength distribution in Legendre polynomials 
directly from the starting vector and Hamiltonian. Wc 
find this method to be stable and effectively reproduce 
distributions to great accuracy. However, the other meth- 
ods discussed here converge more rapidly near the ends of 
the spectrum than in the middle, while the expansion in 
Legendre coefficients does not, in general. The expansion 
also lacks positive-definitcness. 

In the next section, we discuss a fourth method that 
does not attempt to exactly preserve moments, but 
proves in fact to be nearly exact (errors less than 0.01%, 
typically). It is stable, positive definite, and easy to im- 
plement. It is the Lanczos method we recommend for 
those needing to generate response surfaces. 
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IV. THE PIECEWISE MOMENTS METHOD 

The Piecewise Moments Method (PMM) is based 
on the Lanczos matrices and is accurate and positive- 
definite, requiring m + 1 Lanczos calculations to define 
the response function over the entire plane. The basic 
idea behind the PMM is to solve the Lanczos problem 
separately for each of the components \v®), \v{), • • • , \v™) 
of the vector \v±(y)), while incorporating the resolution 
function directly into the algorithm. The resolution func- 
tion provides a prescription for handling inner products 
between Lanczos vectors generated from different start- 
ing vectors. 

While the algorithm is general, it is most transparent 
if formulated first for a Lorentzian resolution function, 
making use of the Lanczos Green's function expansion. 
Combining Eg s . 151 andl 1 71 yields the Lorentzian-smoothed 
response function 



S L {u;,q) = -y 
ir 



T - K e- 2y \c(y)\ 2 ^ K^ihMf 



I— 1 x / 



(21) 



But this can be rewritten as 

S L (v, q) = -y J - K e- 2y {H{y)\4> L {v)) (22) 

where 



\H{y)) = 



l 



uj — H + ia 



[c \v° 1 )+---c m y m K)]. (23) 



The symmetric form of {4>L{y)\4>L{y)) guarantees the pos- 
itive definiteness of the response function. 

At this point we replace the exact Green's function 
in Eq. ESI by its Lanczos equivalent, execute the algo- 
rithm for each separate component of \<j>L,(y)) as a start- 
ing vector, and invoke Eq. 1111 This requires generating 
the m+ 1 Lanczos matrices L(n, \v\}), i = 0, . . . ,m, then 
generating the continued fractions from the entries in the 
tridiagonal matrices, 

G(uj + ia)\v\) -> G n (oj + ia)\v\) 
= c,\{lo + ia)\v\) + ■■■+ gi(u> + ia)|<). (24) 

This gives the PMM expression for the response function 
in the case of a Lorentzian resolution function, 



S L (u,q) = -y J - K e- 2 v V c* Cj y^ 

7T * — ' 



g£(u + i(T)g}(u+ia){vi\v 



(25) 



fej=i 



Note that Lorentzian PMM does not require a diagonal- 
ization of the Lanczos matrices, as the recurrence rela- 
tions for the coefficients g\ use directly the elements of 
the Lanczos matrix |l6j . 



By construction this result has the right form when 
to=0, i.e., when the starting vector has a single compo- 
nent. In this limit (Ufcluj p ) = Ski, and, of course, the pro- 
cedure would yield the exact moments. In more complex 
cases, the resolution function plays an important role in 
interpreting the overlap of Lanczos vectors, as we discuss 
below. 

First, however, it is helpful to generalize the PMM for 
other resolution functions. We use the Gaussian of Eq. [21 
as an example, as other cases are similar: 



SG(u,q) = 



1 



J-K-2y 



(taivMaiy)) (26) 



where 



\MV)) - e-^- H fl^[c»\vl) + ■■■ c m y m \v?)}. (27) 

Now one applies the Lanczos algorithm, evaluating 
L(n, \v\)), i = 0, . . . , m. These m + 1 n-dimensional ma- 
trices are then diagonalized, yielding the Lanczos eigen- 
values {E\} and {'ip i ~ }, / = l,...,n. The appropriate 
complete set can be inserted in each term of Eg. 1271 e.g., 



i=i 



i=i 



(28) 



One thus derives the PMM result 
1 



S G (uj,q) 



(j v 2tt 



a 1 K < 2:1 E < c ^ +3 

i,j=0 



x J2 e -[(-*t) 3 +«-*?)*]/^ (1)^(1)^1^) (29) 
k,i=i 

This "general form" of the PMM requires diagonaliza- 
tion of the Lanczos matrices, just as in Eq.0 Any other 
resolution function can be substituted for the Gaussian: 
the prescription is to "take the square root" of the res- 
olution function, letting it act symmetrically left and 
right, thereby preserving the positive-definiteness of the 
response function. 

An exercise helpful in understanding this result is to 
consider the limit 77 — > N. In this limit the Lanczos 
eigenvectors will converge to true eigenvectors of H, re- 
gardless of the starting vectors. The scalar products in 
Eo. l29l then reflect the resulting orthonormality, indepen- 
dent of the indices i and j. Now consider 77 very close, 
but not equal, to N. In this case one expects, for i =/= j, 
to find nearly identical eigenvectors \ip k ) and \tp J k ) with 
nearly equal but still distinct eigenvalues. If a in the 
resolution function of Eq. \7\ is significantly larger than 
the eigenvalue splitting, the differences will not matter: 
the overlap will be evaluated just as if we had contin- 
ued to the limit 77 — > N. The response function will not 
change as additional iterations are done. Conversely, if 
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a is smaller, the response function clearly will continue 
to evolve, as additional iterations are done, until those 
eigenvalues do become degenerate, on a scale defined by 
a. 

In general, when n -C N, the situation will be more 
complicated, with a number of states within some energy 
range overlapping between Lanczos calculations done 
with different starting vectors, i ^ j. But the Lanc- 
zos algorithm does properly capture the strength for each 
starting vector, omitting high-frequency information. By 
introducing a resolution function directly into the al- 
gorithm we make high-frequency information irrelevant, 
provided n is large enough for the desired a. Thus, even 
though the high-frequency information omitted for i =/= j 
may be somewhat different for these two starting vectors, 
nevertheless one would expect the procedure described 
above to converge. This expectation can be tested nu- 
merically. 



V. NUMERICAL TESTS OF THE PIECEWISE 
MOMENTS METHOD 

To do so we have picked the test case of electromag- 
netic form factors for 28 Si, evaluated in a full sd-shell 
calculation using the Brown- Wildenthal interaction 
The operator notation follows in part Ref. [T^ . 

The Coulomb response function is defined as 



Cj(w,q) 



f 



\{f ; J\\J2Mj(q^)(^^)\\g.s.-,0)\ 2 



(30) 



where 1 1 denotes a reduced matrix element and the sum 
extends over a complete set of sd-shell final states, \ f;J), 
of the requisite angular momentum J. The operator M^ 1 
is 



Mf(qr)=jj(qr)Yj M (n r ) 



(31) 



Because the ground state has J = 0, it is particularly 
simple to rewrite this in the form of Eq. [S] 

Cj(w,q)=J2^- E f)x 
f 

A 



\(f;JM = 0\[J}J2Mj (qn)( 



l + 73(i) 



)\g.s -00)\ 2 



f 



5{uj - E f )\(f; JM = 0|O( 9 )| S . S .;00)| 2 (32) 



where [J] — \J2J + 1. The operator O(q), introduced 
to make the analogy with Eq. [5] clear, can be evaluated 
using the tables of Donnelly and Haxton |l5j . Its second- 
quantized form is 



Mi 

v/47T 



y{J -2 )/2e -yJ2(-iy 

a/3 



jot J 3/3 

-m a mp 



(33) 



where Vm, (y) is the polynomial tabulated in the ta- 
bles. The sums over single-particle states a and [3 are 
restricted to protons because of the isospin projection 
(l+r 3 (i))/2 and to the sd-shell, due to the nuclear model. 
In the case of CO, the inert core nucleons must also be 
included. As the maximum single-particle angular mo- 
mentum in the sci-shell is j=5/2, J=0,2, and 4 Coulomb 
multipoles are possible. 

Similar operators can be obtained for the transverse 
electric response function Ej(uj,q) 



[A 



a 6 



3c 



J 3(3 
lTlf3 



{pf'p) + YPEj(y)) at «p a /3 P + Y p sj(y) at 



(34) 



and transverse magnetic response function Mj(w, q) 



AttM 



y (j-i)/2 e -vj2(-iy 



a/3 



3a 



J 30 
mp 



(vfj{y) - Y^(y)) at "P a /3p - ^-p^vW "„ a /3„J • (35) 

In these equations /i p and \x n are the proton and neutron 
magnetic moments. Normally the charge and magnetic 
single-nucleon couplings are described by form factors, 
but we will treat all couplings as fixed (point nucleon 
limit), as the purpose of this study is the modeling of 
the nuclear momentum dependence. If these couplings 
are given a common momentum dependence, that factor 
could be added to the nuclear results we present below. 

In these equations the underlying single-particle oper- 
ators are 



a/3 , 



(qr) 



Mfj{ q r) ■ iy 
1 



-V x Mfj{qr) 

q 



l - 
-v 
q 



Xf(qr) = Mfj(qr) 
'1 



q 



V x Mfj{qr) 



(36) 



where the spherical Bessel vector harmonic operator is 



Mf L (qr) = 3L(qr)Yfl 1 (n r ) 



M 



(37) 



The polynomials p a ^ can again be found in the ta- 
bles. The overall polynomial behavior - that is, the y- 
dependence other than the overall multiplicative factors 
shown explicitly in Eas. I33I34I35I - are as follows 



CO 
C2/E2 
C4/E4 
Ml 
M3 
M5 



(y\y 2 ,y 3 ) 

(y\y 2 ) 

y 1 

(y°,y\y 2 ) 

(«V) 

y° 



(38) 
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The most complex cases, the CO and Ml response func- 
tions, have three contributing terms in y. Thus a maxi- 
mum of three Lanczos calculations is needed to define any 
electromagnetic response function over the (u>,q) plane, 
for PMM calculations in the sd-shell. 

All of these response functions were evaluated with the 
PMM, for several resolution- function widths (er=1.0, 0.5, 
and 0.25 MeV) and for n ranging from 50 to 400 itera- 
tions. The results we show all assume a Lorentzian for the 
resolution function. The accuracy of the PMM was first 
tested by examining cuts corresponding to y=constant in 
the (u>, q) plane. For such trajectories, 0(y) is fixed, so 
that Eq. [S] can be used - an exact moments treatment. 
This will allow us to test our qualitative argument that 
the PMM should be numerically difficult to distinguish 
from an exact moments treatment, provided n is large 
enough for the given a. Figures 1-3 show the results 
for the CO, C2, E2, Ml, and M3 response functions - 
the cases with more than one component in the vector 
\vi(y)) - for y=2 and a=0.5 MeV, and with n=50,100, 
and 200. 

In those cases where there is significant structure in the 
response function, discrepancies are apparent at n=50, 
but these disappear as more iterations are added. In 
all cases, by rt=200 differences between an exact mo- 
ments calculation and the PMM result are not readily 
discernible on the scale of the graphs. Note that the 
differences at n=50 reflect that fact that neither the ex- 
act moments calculations nor the PMM results are fully 
converged. Similar results were obtained for y=l. This 
kind of test could be made by anyone using the PMM, 
to guarantee, for the chosen resolution function and tr, 
that a sufficient number of iterations have been done to 
produce the desired accuracy. 

The next series of figures provide a more detailed 
look at the convergence and its dependence on a. Fig. 
4 provides benchmarks for exact moments calculations, 
the residual discrepancies between calculations for n=50, 
100, 200, and 400 iterations and one with n=600, which 
we will take as a fully converged result. The figures show 
the convergence in the case of the Ml response evaluated 
along the line y=l at the resolutions a = 0.25 and 1.0 
MeV. The behavior is just as one would expect. The 
missing contributions are oscillatory, with a frequency 
that is roughly proportional to n but virtually indepen- 
dent of a - , for the range we explored. The envelope of 
the oscillations decreases with increasing n, shrinking by 
more than an order of magnitude for every additional 50 
iterations for cr=1.0 MeV. The decrease slows to about 
a factor of 5 every 100 iterations for the more taxing 
calculation with cr=0.25 MeV. 

Similar calculations, not shown, were performed for the 
CO response, which has much less structure than the Ml 
response. The results are qualitatively similar to that 
shown in Fig. 4, though the convergence of the envelope 
is a factor of ~ 30 more rapid. 

Figure 5, the analog of Fig. 4, gives the Ml PMM 
residuals (again relative to an exact moments calculation 
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FIG. 1: Comparison of an exact moments calculation (solid 
line) with the PMM (dashed) method for the CO and C2 re- 
sponse functions along the y=2 line in the response plane. 
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FIG. 2: As in Fig. 1, only for the E2 response. 



with n=600). While again an oscillatory pattern emerges 
with a frequency like that of Fig. 4, its structure is less 
regular. This is the result of the interference between 
the three Lanczos patterns that contribute to the PMM 
result, corresponding to the starting vectors \v\), 
and \vf). The PMM envelopes tend to be a factor of 
two to three more extended than those from the exact 
moments calculation, though in one case the difference is 
larger. 

Figure 6, the residuals for the PMM calculation for 
the CO response, with y=l and a = 1.0 MeV, show an 
interesting effect. Through most of the spectrum the 
same oscillatory features and diminishing envelope with 
increasing n are seen. But in the low-energy region, a 
persistent feature has converged by n=100. It appears 
dominantly positive, but as the fractional convergence is 
graphed and as this response is dominated by the ground 
state, where the differential is negative, this is mislead- 
ing. The low-energy CO response is characterized by iso- 
lated extremum eigenvalues, as is apparent from Fig. 1, 
with the ground state carrying the entire response from 
\ v i) (v ~~ * 0)- These appear to be conditions that al- 
low for a small deviation of the PMM from the results of 
an exact moments calculation, in converged calculations. 
However, the deviation is small, less than 0.01%. Very 
similar effects were found in the CO responses for a = 
0.5 MeV and 0.25 MeV, with the discrepancies reaching 
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FIG. 3: As in Fig. 1, only for the ml and M3 responses. 
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FIG. 4: Comparison of the convergence of the Ml response FIG. 5: As in Fig. 4, only for the PMM calculation, 

for j/=l for a — 0.25 and 1.0 MeV in an exact moments cal- 
culation. 
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FIG. 6: Convergence of the PMM CO response for y=l and 
cr=1.0 MeV. A small difference persists near the dominant 
extremum eigenvalues. 



0.03% and 0.06% in these cases, respectively. 

Finally, we show a series of PMM results for the re- 
sponse surface. Figure 7 shows contour plots for the CO, 
C2, and C4 responses, while Figs. 8 and 9 give the three- 
dimensional projections of the E2 and E4 and the Ml, 
M3, and M5 responses, respectively. The general shift of 
strength to larger y with increasing multipolarity is ap- 
parent. In most applications of the PMM, such response 
surfaces, determined as a function not only of (ui,q) but 
also of the oscillation parameter 6, would be the end re- 
sult of the calculations. 



VI. CONCLUSIONS 

An attractive property of the matrix elements of elec- 
troweak single-particle operators between harmonic os- 
cillator states is that they can be evaluated analytically, 
yielding a simple form that includes a finite polynomial in 
y = (qb/2) 2 . This property has been exploited frequently 
in calculations of discrete transition amplitudes. In this 
paper we have shown that moments methods exploiting 
this property can very efficiently characterize the entire 
response surface over the (u>, q) plane. This can be viewed 
as an important extension of well-known moments tech- 




Si C2: 200 iterations, <r L =0.5 MeV 




Si C4: 200 iterations, <r L =0.5 MeV 



g> 15 
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FIG. 7: The PMM response surfaces for the CO, C2, and C4 
multipoles, with contours drawn at successive factors of 0.5 
of the maximum, until 0.001 is reached. 
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FIG. 8: The PMM response surface for the E2 and E4 multi- 
poles. 



niques for determining the distribution of Gamow- Teller 
strength along the line q = in the (to, q) plane. 

While the method was motivated by the observation 
that moments must have a polynomial form in y, algo- 
rithms we designed to exactly preserve the moments were 
found to be numerically unstable, if a sufficient number 
of iterations n were done. This difficulty is the well- 
known "classical moments problem," the determination 
of a discrete distribution from knowledge of the distribu- 
tion's moments. An alternative method, called the Piece- 
wise Moments Method, was introduced in which each 
polynomial component of the starting vector was treated 
separately in the Lanczos procedure, while also incorpo- 
rating a resolution function directly into the algorithm. 




FIG. 9: As in Fig. 8, only for the Ml, M3, and M5 responses. 
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The method is extremely stable, positive definite, and 
effectively equivalent, numerically, to an exact moments 
reconstruction. After convergence - which typically re- 
quires ~ 200 iterations for a resolution a ~ 0.5 MeV - we 
found maximal differences between the PMM and exact 
moments calculations of about 0.01%, in one case (the 
CO multipole). The number of iterations n required for 
convergence increases with decreasing a. 

We noted that the absence of a moments method to re- 
construct the response in the full (q, w) plane had previ- 
ously led to the use of very simple nuclear models, so that 
state-by-state sums over transitions could be performed. 
Clearly the PMM will now allow theorists to perform 
analogous calculations using state-of-the-art shell-model 
wave functions in very large model spaces. The efficiency 
with which the PMM constructs the response function 
over the response plane is impressive. In the example we 
explored - the electromagnetic form factors for 28 Si in 
the sd-shell - the most complicated multipoles, CO and 
Ml, required only three Lanczos steps. 

Another aspect of the method, discussed in Section II, 
is that it can be viewed as a numerical effective theory 
in the sense that exactly that information needed to re- 
construct S(u>,q) can be systematically extracted from 
exceedingly complicated nuclear structure calculations. 
Specifically, if n is the number of iterations and m the 
rank ofp(y), one needs (m+ l)(2n— 1) tridiagonal Lanc- 
zos matrix elements and [m(m — l)/2]n 2 Lanczos vec- 
tor overlaps to implement Eq. 1251 for example. Once 
this information is in hand, a simple routine could be 
coded to generate S(u>, q) as a continuous function of ui 
and q, as well as of the oscillator parameter b and a. 
(The one caveat in that for any fixed n, there will be 
some minimum a beyond which the number of iterations 
n performed would not be sufficient.) In effect, S(u>,q) 
would be no more complex, in numerical calculations, 
than some analytically known leptonic scattering kernel. 
One application we have in mind is neutrino reactions in 
core-collapse supernova modeling. This approach would 
allow one to model such reactions with state-of-the-art 
shell model techniques, yet produce a result sufficiently 
simple that it could be used on-line within a sophisti- 
cated supernova code. The fact that S(ui,q) is given as 
a function of u and q, rather than as a grid of values, 
is also important in such applications. This will allow 
supernova modelers to more easily guarantee properties 
such as detailed balance in reactions and inverse reac- 
tions - which if not enforced exactly, can lead to energy 
generation and other spurious physics when weak inter- 
actions are in equilibrium. Detailed balance can be more 
difficult to enforce in cases where S(u>,q) is provided on 
some discrete numerical grid. 

We would like to mention three follow-up studies that 
we think will make the present work more valuable. First 
is the extension to weak interactions. This is relatively 
simple, as only two new operators (in addition to those 
we have treated in Eas. l31l and l35|) arise when axial cur- 
rents are added (in the standard nonrelativistic treat- 



ment where charges and currents are kept to order 1/M) 
[LH fl*5| . One of the main motivations of the present 
paper is to develop a technique that can be applied to 
neutrino response functions in the energy range up to ~ 
1 GeV. 

A second question is the loss of unitarity for 
momentum-dependent operators acting in finite shell 
model spaces. This issue does not arise for the standard 
application to the Gamow- Teller operator g e /^ &t±, for 
any complete nhui shell-model space, as the operator has 
no cross-shell matrix elements. But for more general op- 
erators, as the momentum increases, so does the strength 
of excitations to states outside the model space. For a 
vector 0(y)\g.s.) it is simple exercise to calculate the loss 
of probability due to excitations outside the model space. 
If one fails to take into account such loss of probability, 
trends in inclusive cross sections as a function of q will be 
distorted. Thus it is important to either correct for such 
effects, or to evaluate their size to estimate uncertainties 
in results. 

A third issue is the overcompleteness of shell-model 
spaces due to spurious center-of-mass motion. This can 
be troublesome when dealing with one-body operators 
that can excite center-of-mass excitations. While this is- 
sue did not arise in our example of 28 Si because we con- 
fined ourselves to the hu> sd-shell space, it will in more 
complicated spaces. If the space is separable - e.g., any 
nhio shell-model calculation with oscillator wave func- 
tions - the standard technique, when dealing with dis- 
crete calculations, is to remove spurious states by adding 
to the Hamiltonian a term aHcM, where Hqm is the cen- 
ter of mass Hamiltonian and a ~ 100 a coefficient cho- 
sen to "blow out" spurious states from the low-energy 
spectrum [Tlj . This works well for converged Lanczos 
calculations: the addition of such a term forces low- lying 
eigenvalues to have the center-of-mass in the Is state. 
One would need to explore numerically whether a simi- 
lar technique might allow some approximate separation 
of spurious excitations over the full spectrum: we are not 
aware of any studies of this method apart from the case 
of converged extremal eigenvalues. Alternatively, if the 
shell-model Hamiltonian is (properly) translationally in- 
variant and \g.s.) has been constructed so that its center 
of mass is in the Is state, center-of-mass excitations in 
the vector 0(y)\g.s.) could be removed at the outset. 
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APPENDIX: ALTERNATIVE METHODS 

In this Appendix we discuss in more detail alterna- 
tive Lanczos response function methods that we have ex- 
plored. This discussion might be useful to those that 
would like to further explore some of the stability issues 
we encountered. The first two approaches have as their 
goal a reconstruction of S(lj, q) that exactly captures the 
information in the 2n — 1 moments. The main drawback 
in both methods is instability for moderate n, connected 
with the classical moments problem (the inversion from 
moments to a distribution). A third method is discussed 
which deals directly with distributions while also preserv- 
ing a specified number of moments. 

1. Naive Moments Method 



While sound mathematically, the NMM is problem- 
atic numerically because the solution to the classical mo- 
ments problem embodied in Eqs. IA.2I and I A. 31 is not a 
stable one. For example, Whitehead and Watt provide in 
Ref. |l<Sl a simple 4x4 matrix example in which signif- 
icant loss of accuracy occurs. There have been attempts 
[Tfll |23| to finite alternate methods that use moments in 
a manner that will improve the inversion to a distribu- 
tion. None of these have proven to have the stability of a 
direct Lanczos construction, however. In the numerical 
tests we performed, significant errors typically occurred 
for n ~ 20, in calculations performed with 64- bit accu- 
racy. As is apparent from many of the calculations pre- 
sented in this paper, often 50-200 iterations are required 
before a response function is accurately reconstructed, 
for resolutions we explored. 



As discussed in the main body of the paper, the NMM 
begins with a starting vector of the form 

c{y)\ Vl (y)) = coK) + c iy \v\) + ... + c m y m \v™) (A.l) 

from which the moments (vi(y)\H x \vi(y)), A = 
1, . . . , 2n— 1, can be determined for any y once the mixed 
moments have been evaluated. As knowledge of the 
moments as a function of y is mathematically equivalent 
to knowledge of L(n, \ v±(y))), the response function can 
thus be evaluated as a function of both y and to. 

On first sight, it appears that the inversion from mo- 
ments to the tridiagonal matrix can be easily done. It has 
been shown ^3] that the elements of the Lanczos matrix 
are related to the moments by determinants so that 



and 



= Mi-i/£i-i -M-2/A-2 (A.2) 



ft = C\ /2 C]^/ C t -i, (A.3) 



where the determinants C and M are defined by 



1 Ml M2 
Ml M2 M3 



and 
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1 Mi M2 
Mi M2 M3 

Mn • • • 



Mn 
Mri+1 

M2n 



Mn-1 Mn+1 
Mn Mn+2 



M2n-1 M2n+1 



(A.4) 



(A.5) 



with fix = ( v i\H x \vi), and with A / f_i=0, Mo = Mi, an d 
Cq=1. The NMM uses these equations to determine the 
Lanczos matrix L(n, \vi(y)} from the fix{y)- Then one 
can proceed in the usual way to diagonalize this matrix 
and then construct the strength function, using Eq. 



2. Legendre Moments Method 

Operationally, the rapid loss of precision in the NMM 
occurs because of severe cancellations in the determinants 
fEas. lA.4l and lA~5)l that arise from the dominance of the 
largest eigenvalues in the highest-order moments. This 
dominance occurs because naive moments methods work 
with a set of non-orthogonal basis vectors: H n \vi) with 
n non-negative integers. The Lanczos algorithm, on the 
other hand, builds a set of orthogonal vectors \v n ) as it 
goes, each carrying information about H that is indepen- 
dent from that in its predecessors, and this produces the 
great stability of the algorithm. 

We attempted to find a new way to iterate for ai and 
ft in terms of moments that, like the original Lanczos 
algorithm, are built on a set of orthogonal vectors, \u)i), 
constructed as part of the iteration. We describe a first 
attempt at an iterative method, which fails because it is 
equivalent to the NMM. We then give a closely analogous 
procedure based on linear combinations of moments that 
proved, at times, to be considerably more stable. 

We begin as in the NMM procedure by calculating the 
mixed moments O^-, from which we can then evaluate 
the moments (vi(y)\H m \vi(y)} as a function of y. 

We then build up the orthogonal vectors \wi) of dimen- 
sion i iteratively, starting with 



Wi 



W 2 



(«i) 



1 

ft 



(A.6) 



(A.7) 



where a\, ft are easily identified from the moments 
(vi\H\vi) = a.\ and (vi\H 2 \v 2 ) — a\ = fi\ (defined as 
in ordinary Lanczos). With these vectors in hand we 
compute the next on, ft, and \w n ) from 



— 2(ai + a 2 ■ 



+ a r , 



(A. 
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where L n is the truncated Lanczos matrix of Eq. [3 

In this way we construct the Lanczos matrix, which 
can then be diagonalized and used to derive the response 
function in the usual way. This algorithm yields signifi- 
cant stability improvements over the procedure outlined 
in the NMM discussion. However, the new algorithm still 
suffers from cancellations between terms in both Eas. IA.8l 
and IA.9I Successively higher moments yield increasingly 
large numbers on each side of the minus signs in those 
equations. These numbers are subtracted to yield suc- 
cessive ati and (3i that remain of order unity. Such large 
cancellations lead to loss of precision that, though less se- 
vere than in naive moments methods, has the same root 
and the same consequence of failure in that eventually 
< 0. At this point the algorithm fails. 
One possible solution to this problem is to find combi- 
nations of moments that remain relatively stable in size, 
e.g., some set of orthogonal polynomials. If we scale 
and shift energies so that the range of eigenvalues can 
be mapped onto [-1,1], one obvious choice would be Leg- 
endre polynomials. They contain the same information 
as the moments, and their matrix elements will have the 
same dependence in y. However, as they have a fixed 
magnitude at the endpoints, they are not overwhelmed 
by the extremal eigenvalues for large n. In the LMM the 
moments of the NMM are replaced by 



(vi\Pl\vi) 
(v 1 \P 2 HP 2 \v 1 ) 



{qaw) 

(qi\H\qi) 
(92192) 

(q2\H\q 2 ) 



The operators Pi (H) are linear combinations of powers of 
H, with coefficients identical to those on the correspond- 
ing powers of the scalar x in the Legendre polynomi- 
als Pi(x). Instead of computing moments, we calculate 
the overlaps (qi\qi), using the Legendre polynomial re- 
currence relation to generate \qi). The vectors \qi) have a 
simple expansion in y because they are built from linear 
combinations of the \v]), so 



l«) = \qf: 



\qf)y 



\Ql)V 



(A.ll) 



and \q\) is easily computed from \v\) at the start of the 



calculation. The recurrence relations used are then 
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It is straighforward to find the Lanczos matrix L n in 
terms of the (q n \q n }. This resulting inversion - from Leg- 
endre moments to the tridiagonal matrix - proved to be 
significantly more stable than the NMM inversion. The 
recurrence relations to compute L n are 



{q n -x\H\q n _\) - (w n -x\L n -x\w n -i) 



(2n-3)!! 
("-1)1 



(a<i +Q< 2 ■■■ + £*„) (A.15) 
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(A.17) 



with \wi), \w 2 ), Q.1, and pi as before. 

These two moments methods (Eqs. IA.8HA.10l and 
Eqs. IA.15llA~T7|) are defined by similar procedures. We 
compare their stability in Fig. 10. In our experience, the 
Legendre moments algorithm can sometimes run to very 
high numbers (> 80) of iterations with no significant loss 
of precision, but it also sometimes runs aground quickly 
on cancellations in Eq. IA. 161 The lack of predictability 
is clearly an issue. 

Our failure to identify a more reliable method for deter- 
mining distributions from moments motivated us to look 
for other approaches, e.g., ones that might not capture 
exact information on the moments at every y or exactly 
guarantee positive-definiteness of the response function, 
but would remain stable and accurate under continued 
iteration. This led to the PMM method we favor, as 
well as one other moment-preserving approach described 
below. 



3. Legendre Coefficients Method 

Rather than constructing a Lanczos matrix, diagonal- 
izing, and finding the strength distribution, an approx- 
imate strength distribution can be computed directly 
from an expansion in Legendre polynomials, 

00 

S(cj, y) = y J ~ K e- 2 y £ ai {y)Pi{u). (A.18) 
1=0 

The coefficients ai (y) can be computed using the orthogo- 
nality relation for Legendre polynomials and the identity 



y J - K e- 2 y\c{y)\\v 1 \P l {H)\v 1 ) = I S(u,y)P t (uj)dcj 



(A.19) 
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FIG. 10: The number of iterations before loss of precision 
produces /3i < in the moments algorithms, for each operator 
computed in our 28 Si example functions of momentum 

transfer. The upper graph is the result for the method of 
Eqs. IA.6HA.10l , while the bottom graph gives the Legendre 
Moments Method fEas. IA.15HA.17J . 



so that 
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We refer to the method of computing response functions 
from Legendre coefficients computed with these relations 
as the "Legendre coefficients" method. 

Here again, H must be shifted and scaled to the inter- 
val [—1,1] before calculating these matrix elements, then 
shifted and scaled back at the end to obtain the true 
strength distribution. A numerically stable approach to 



calculating the ai{y) uses the same vectors \q n ) defined 
for the Legendre moments method, constructed in the 
same way, so that 



ai(y) 



21 + 1 



j,fc=0 



CjC k y [v 
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where Cj are defined in Eq. IA.1I 

The approximate strength function S^c represented 
by the first n Legendre polynomials will always contain 
oscillations about the true strength function, at an en- 
ergy scale given by the order of the highest Legendre 
polynomial in the expansion, so roughly on the scale 
Aw/n, where Auj is the difference between the maxi- 
mum and minimum eigenvalues of H. (Since the true 
response function is a sum of delta functions, its expan- 
sion in Legendre polynomials never truncates.) There 
is no restriction on the sign of S^ci^, y), so its oscilla- 
tions may make it negative in places even though S(u), y) 
is strictly positive. For these reasons, practical use of 
this method would probably require smoothing smooth- 
ing the function S^ci^i v) 1,0 finite resolution over scales 
of approximately Atu/n. Since we have a Legendre ex- 
pansion of S 7 l Cl the convolution to produce the smoothed 
function may be performed efficiently (i.e., reduced to a 
series of matrix multiplications) by working with Legen- 
dre expansions of the smoothing function. The Legendre 
coefficient expansion may be useful in applications where 
the response function needs to be convolved with some 
other function, for example, a thermal energy distribu- 
tion. On the other hand, while the Lanczos methods 
converge more rapidly at the ends of the eigenvalue spec- 
trum than in the middle, that is not in general the case 
for the Legendre expansion. 

We also note that, whereas the ordinary Lanczos al- 
gorithm and our variants on it need to run n iterations 
on each vector piece to reproduce 2n — 1 moments of 
the response function, the Legendre coefficients methods 
needs to run 2n — 1 iterations to produce 2n — 1 mo- 
ments, essentially twice as long as the other methods in 
this Appendix. As in the Lanczos method, each iteration 
contains as its basic time-consuming step a matrix- vector 
multiplication in the original large basis. 

It has come to our attention that a closely-related 
method is in use in physical chemistry, where it is used 
to compute the quantum-mechanical time evolution of 
molecular states [21], |22j . Authors working in that area 
point out that if Chebyshev polynomials are used instead 
of Legendre polynomials, the vectors \qi) may still be 
found by recurrence, but only the first n such vectors 
need to be computed in order to find the first 2n terms 
of the expansion in Chebyshev polynomials. 
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